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Abstract 

We give a detailed comparison of the hierarchical quantum master equation (HQME) method 
to a continuous-time quantum Monte Carlo (CT-QMC) approach, assessing the usability of these 
numerically exact schemes as impurity solvers in practical nonequilibrium calculations. We review 
the main characteristics of the methods and discuss the scaling of the associated numerical effort. 
We substantiate our discussion with explicit numerical results for the nonequilibrium transport 
properties of a single-site Anderson impurity. The numerical effort of the HQME scheme scales 
linearly with the simulation time but increases (at worst exponentially) with decreasing temper¬ 
ature. In contrast, CT-QMC is less restricted by temperature at short times, but in general the 
cost of going to longer times is also exponential. After establishing the numerical exactness of the 
HQME scheme, we use it to elucidate the influence of different ways to induce transport through 
the impurity on the initial dynamics, discuss the phenomenon of coherent current oscillations, 
known as current ringing, and explain the non-monotonic temperature dependence of the steady- 
state magnetization as a result of competing broadening effects. We also elucidate the pronounced 
non-linear magnetization dynamics, which appears on intermediate time scales in the presence of 
an asymmetric coupling to the electrodes. 

PACS numbers: 85.35.-p, 73.63.-b, 73.40.Gk 
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I. INTRODUCTION 


Impurity problems are ubiquitous in the theoretical description of nonequilibrium systems 
BE3- They constitute small entities with a limited number of degrees of freedom that are 
coupled to reservoirs with continua of non-interacting degrees of freedom. One intuitive 
physical realization of such a model is a molecule adsorbed on a surface or contacted by 
electrodes [3]. A variety of nonequilibrium scenarios may be described in terms of impurity 
models, for instance by preparing the impurity in an excited initial state or by coupling 
it to different reservoirs in different thermodynamic states. Another intriguing application 
occurs in dynamical mean field theory [2| 0H6j , where lattice problems either in or out of 
equilibrium are mapped to impurity problems with an environment that is determined by 
a self-consistency criterion. This has been important, for instance, in understanding the 
metal-insulator transition in materials like transition metal oxides HIM and has become 
an important paradigm in studying nonequilibrium effects in extended interacting systems, 
including thermalization after an interaction quench mi, the nonequilibrium steady state 
mm and Bloch oscillations [5j [12], 13] under the influence of a static electric field. Thus, 
the theoretical description of impurity problems is a key element in understanding a wide 
range of phenomena, in particular nonequilibrium effects. 

Few exact solutions are available, and a number of methods have been developed in the 
past decades to solve nonequilibrium impurity problems. They can be sorted into two broad 
categories: approximate and numerically exact methods. Typically, numerically exact meth¬ 
ods allow us to simulate some property of a model in what might be considered a numerical 
experiment. Approximate methods, on the other hand, may miss important physics or suffer 
from artifacts due to the approximations involved. A combination of methods, which oper¬ 
ate on different levels of approximation is often useful and helps to elucidate the relevant 
physical mechanisms jTT EU] . 

The nonequilibrium Anderson impurity model has been treated by several numerically 
exact methodologies. Some approaches require a discretization of the electrodes, for example 
density-matrix renormalization group |2Tl - l26] . numerical renormalization group pH EH EE] 
or multilayer multiconfiguration time-dependent Hartree theory [2HIEUJ. These methods are 
useful at low temperatures and/or voltages. They are restricted by revival oscillations and 
a limited spectral resolution of the leads [31]. Other approaches can take advantage of the 
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noninteracting nature of the leads and do not require discretization. This includes iterative 
path-integral schemes |32HS^, which converge only for a limited set of parameters, and 
stochastic schemes [1171 ITT] , where the growth of the statistical error restricts the accessible 
time scales. In the presence of a short memory timescale, long time scales can be accessed 
by a combination of reduced dynamics techniques with a short-time numerically exact 
scheme; or by the hierarchical master equation method m where the numerical effort 
scales linearly with the simulation time. The latter, however, can only be converged if the 
temperature in the electrodes is not too low mm- 

The numerical effort associated with most numerically exact schemes restricts practical 
calculations to specific limits |36] or limited ranges of parameters. It is important to delineate 
the regimes in parameter space to which each method is applicable, and in particular to find 
out where exact results are not available. In this work, we elucidate the practical limitations 
of two numerically exact schemes: the continuous-time quantum Monte Carlo (CT-QMC) 
method [3T;, 39], HU, 13] 0H H3-09] and the hierarchical quantum master equation (HQME) 
method HOT EM. We will discuss the main features of these approaches, characterizing, 
in particular, the associated numerical effort. We find that the range of parameters where 
the two methods can be applied overlaps, but also exhibits areas where only one of the 
methods can be applied. As we will see, HQME turns out to be the method of choice to 
study long-term dynamics if the temperature of the reservoirs is not too low. In contrast, 
CT-QMC gives access to the short- and intermediate-time dynamics over a wider range of 
temperatures. 

We demonstrate our findings using an archetypal nonequilibrium problem: transport 
through a single-site Anderson impurity that is coupled to left and right electrodes (see Fig. 
[I] for a graphical representation). The most obvious physical realizations of this impurity 
problem are quantum dots containing a single spin-degenerate level. The first such real¬ 
izations were based on quantum-confinements in patterned semiconductor heterostructures 
[53) El]. Single-molecule junctions often exhibit similar behavior [55H58] . but in a setting 
where experimental techniques give less control over the parameters of the junction. Addi¬ 
tionally, other effects, e.g. due to vibrational degrees of freedom, are important [59H69] . In 
these systems, transport is induced by shifting the electrochemical potentials in the leads 
with respect to each other such that electrons tunnel through the impurity in order to move 
from the lead with higher chemical potential to that with lower chemical potential. In semi- 
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conductor heterostructures this shift is achieved by charging or discharging the leads (i.e. 
by filling or emptying electronic levels; cf. Fig. [TJd) . In single-molecule junctions, the leads 
are less likely to be charged, and the shift of the electrochemical potentials is accompanied 
by a shift of the respective conduction bands (cf. Fig. []J). We will show that these different 
ways of inducing transport strongly affect the initial dynamics of the impurity. 

We also present exact results for the complex magnetization dynamics of an Anderson 
impurity in various nonequilibrium situations. Typically, this quantity exhibits the slowest 
relaxation behavior (43J and, as we will see, exhibits a non-linear behavior on all times 
scales, in particular when an asymmetric coupling to the electrodes is considered. To date, 
this dynamics had only been accessible at great computational cost using state-to-the-art 
CT-QMC methods combined with reduced dynamics [43]. The HQME method gives access 
to exact results of this long-lived correlated dynamics and allows us to perfrom a scan over a 
wide range of parameters. It can also be used to derive approximate results and, therefore, to 
study the influence of higher order processes. We are, therefore, able to elucidate the origin of 
the non-monotonic temperature dependence of the magnetization that was recently reported 
in Ref. [03] to be the result of competing broadening effects. We also End a pronounced 
non-linear behavior of the magnetization on intermediate (still rather long) time scales (cf. 
e.g., Figs. [7] and [9]). In passing we note that the nonequilibrium Anderson impurity model 
and its generalizations are of great interest in the field of strongly correlated materials within 
the dynamical mean held theory approximation la m\. 

The outline of the article is as follows: I 11 Sec. [TTJ we present the theoretical methodology. 


This includes a short description of the single-site Anderson impurity model (Sec. IIA), 


the HQME method (Sec. IIB) and the CT-QMC approach (Sec. IIC). A discussion of 


practical aspects of the two methods is given in Sec. IIP Numerical results on the time- 


dependent transport properties of an Anderson impurity are presented in Sec. Ill where 
we first formulate the different ways of inducing transport and detail the model parameters 
(Sec. |III|A). A direct comparison of results that are obtained by the HQME method and 


the CT-QMC approach is presented in Sec. Ill B, where we follow the time evolution of the 
electrical current that is flowing through the impurity, starting from a product initial state 
where the impurity is not populated by electrons. These results represent the first explicit 
validation that the HQME approach gives numerically exact results. We also explore how 
the choice of whether or not to shift the conduction band with the applied bias voltage affects 
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FIG. 1. Panel (a): Graphical representation of an Anderson impurity, which is realized by a 
quantum dot (QD). The dot is coupled to a left (L) and a right electrode (R). Panel (b): Single¬ 
particle levels of the quantum dot junction depicted in Panel (a). The shaded areas depict the 
occupied states in the electrodes, which, for simplicity, are assumed to provide Lorentzian-shaped 
conduction bands. Applying a bias voltage to such a system means that the corresponding chemical 
potentials /xl/r are shifted and the electrodes become charged. Panel (c): Single-particle levels of a 
molecular junction. In contrast to the quantum dot realization (Panels (a) and (b)), the conduction 
bands are shifted in the same way as the chemical potentials. Thus, applying a bias voltage, the 
electrodes do not become charged. 


the results. We then discuss the magnetization dynamics of the impurity in the presence 
of an external magnetic held, considering both a symmetric and an asymmetric coupling to 


the electrodes (Sec. Ill C). 


II. THEORY 

A. Model Hamiltonian 

We study the transport properties of an Anderson impurity that is coupled to a left (L) 
and a right (R) electrode or lead. The Hamiltonian of this well established system, 


H — //imp + //[. + Hr + //turn 


( 1 ) 
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can be decomposed into the impurity Hamiltonian, fA mp ; the left and the right lead Hamil¬ 
tonians, and II\\ : and a coupling operator H t un . The impurity Hamiltonian 

I E e a d\d a + Udjdfdjdj.. (2) 

<re{t,4-} 

represents an electronic level that is addressed by creation and annihilation operators d\ 
and d a . It can hold a single spin-up (d) or spin-down (d) electron at energies and e^, 
respectively, where = ej. without the influence of an external magnetic held. It can also 
hold a spin-up and a spin-down electron simultaneously. Such double occupation costs an 
additional charging energy U > 0, representing repulsive Coulomb interactions. 

Each lead is described by a continuum of non-interacting electronic levels 

Hh/R — 'y ^ ^akC a f.C a k (3) 

L/R,cr 

with energies e a k- These levels are addressed by annihilation and creation operators c a k 
and c\ k . The coupling between the impurity and the electrodes can be characterized by 
tunneling matrix elements V a k, and the corresponding coupling operator is written as 

H tun = (V*kct k d a + h.c.). (4) 

/cE'fLH-R} ;cr 

The resulting tunneling efficiencies, or level-width functions, 

r*>(e) = 2nY J Kk\ 2 5{e-e ak ), (5) 

k&K 

depend on the energy of the tunneling electrons (K e {L,R}). Throughout this work, we 
assume that the left and the right electrode have the same properties, in particular that 
they have the same temperature T. The only difference between the electrodes occurs in 
the presence of a bias voltage $ = — /ir ^ 0. The position of the chemical potentials 

of the left and the right lead, /xl and n r, are shifted in different directions. We assume a 
symmetric voltage drop such that [M, = <f>/2 and fi r = —<f>/2. The latter is an assumption 
that, however, is not crucial for our discussion. 

B. Hierarchical master equation approach 

We use two different methods to obtain the transport properties of an Anderson impurity. 
The Erst of these methods is the HQME approach [191 201 150H52] . The second one is the 
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CT-QMC method [S3 USE 0310U 08J 09] . For completeness and to establish the notation, we 
outline the basics of the HQME method in this section. The CT-QMC technique is detailed 
in the following section, Sec. 11 C| 

The central quantity of the HQME technique is the density matrix of the impurity 


P = Y1 pl ’ 1 ' iV’imp.i) I (6) 

i,v 

where the V ; im P ,z represent the corresponding Hilbert space. It is obtained by solving the 
hierarchy of equations of motion 




,(«) 


‘impj rji..j 
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J«+l j^Vc+l 


E , ,»A „(«) 
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Ae{i..«} 

VK^*AP^/ j ^ t )+ J2 

AE{ 1 ..«} 

(k+1) 


(t) (7) 

(-1 ) X V’kC,^Ai~S/3x^ 




d SK+1 o- ’■ — 

u, a K+ iPji..] K j K+1 \ L ) \ Pji-JkJk+i 


si) 


A detailed derivation of these equations can be found in Refs. HSH 52 J. The density matrix of 
the impurity enters at the 0 th tier as p^(t) = p(t)- The auxiliary operators Pj*..j K {t) encode 
the dynamics of the impurity that originates from the coupling to the electrodes, starting 
from a product initial state or, equivalently, P^[ Jk (0) = 0 (k > 1). They are distinguished 
by superindices j\ = (K, a, s,p ), which involve a lead index K, a spin index a and an index 
s G {+, —} that corresponds to the creation and annihilation of electrons. The index p is 
related to a decomposition of the lead correlation functions 

den 


C s k,M= I 






( 8 ) 

(9) 


by a set of exponential functions, e _a Ap*, where we use the short-hand notations (cn) = 
/x(cn) and = 1 —//< (cn) with fx{ uj) representing the Fermi distribution function of lead 

K . The use of exponential functions facilitates a systematic closure of the hierarchy (Jt|) [[191 
[52]. We obtain the frequencies c o s Kp and the amplitudes , using a Pade decomposition 
scheme [70], [Zlj • Explicit expressions can be found in Refs. [19] and [20] 

In principle, the solution of the full hierarchy ([7]) is exact. In practical calculations, 
however, the number of Pade poles that can be included is limited. For the present studies, 
we obtained converged results using 100 Pade poles. In addition, the hierarchy of equations 
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of motion (Jt|) needs to be truncated. To this end, we estimate the importance of the operators 
by assigning them the following importance value: 

1 


n 


n 


^K x ,a x ,Px 


( 10 ) 


,A=iL ^A'=1..A Re K,ft']/ \A=1 A ./c . 

We include only those operators in our calculations which have a value larger than a certain 
threshold value A t h- In addition, all operators of the zeroth and first tier are taken into 
account regardless of their assigned importance value. The truncation allows us to reduce 
the numerical effort to a practical level. While this procedure represents an approximation, 
exact results can be obtained by systematically reducing the threshold value A t h until the 


results converge to within the desired precision. As the importance criterion (10) involves 
the ratio between the amplitudes VK,a, P an d the frequencies w s K , that is, effectively the 
ratio Ta' )(T /T, convergence can be achieved more easily at higher temperatures, and the 
numerical effort increases substantially at lower temperatures. We will elaborate on this 


statement in Sec. Ill B , where we show that “large enough” means in the present context 
that the temperature should be above the Kondo temperature. An example of a convergence 
analysis is given in the appendix. 

A central characteristic of the technique is that the equations of motion ([7]) are local in 
time. Thus, the numerical effort of computing the time-dependent operators Pjlj K (t) scales 
linearly with the simulation time t. All of our numerical evidence (cf., e.g ., the appendix) 
shows that the quality of the associated effective expansion of the time evolution operator 
is independent of the simulation time t. We can therefore conclude that the numerical 
effort of the HQME scheme scales linearly with the simulation time. This allows us to 
access the nonequilibrium dynamics of an interacting impurity system on extremely long 
time scales (see, e.g., Ref. ra, where we simulated the dynamics of an interacting quantum 
dot system up to t ~ 10 4 /T). The price of this is a large number of unknown auxiliary 
operators ( t ) that need to be determined. At a given time t , they encode the history 
of the interplay between the impurity and the electrode at earlier times, and contain all the 
information necessary to continue propagating the density matrix to the next time step. 


In addition to the importance criterion (10), the hierarchy can be truncated at a specific 


tier R. This corresponds effectively to a hybridization expansion of the time evolution 
operator of the density matrix p{t) that is valid up to 0{Y nK ’ a /mii\(D,k^T) K ) [T9] , Such 
a truncation is not exact but facilitates a perturbative analysis, which, in principle, can 









be driven to arbitrary order. We can therefore assess the importance of each tier/order by 


comparison to the exact converged results (see Sec. |III C[ ). In this context, it should be noted 

. In the 


that the hierarchy (Jt]) terminates automatically at the second tier for U —>■ 0 
appendix, we demonstrate the convergence of our approach to the exact U = 0 result. 


C. Continuous-time quantum Monte Carlo approach 

In order to establish the numerical exactness of the HQME approach, it is useful to 
compare it to another numerically exact approach based on entirely different principles. 
As noted in the introduction, a wide variety of numerically exact methods with various 
advantages and limitations has been applied to nonequilibrium impurity models [151 IT51 [Ml 
E3 E21 m, HD [73H79]. We have chosen to compare our results with those of a continuous¬ 
time quantum Monte Carlo method [751180]. CT-QMC algorithms are capable of solving a 
variety of impurity models by stochastically summing all terms in an exact diagrammatic 
expansion around some analytically solvable limit. 

Dynamics and nonequilibrium require a real time (rather than an equilibrium imaginary 
time) formulation of the method to conserve numerical exactness. The first real-time imple¬ 
mentations addressed vibrations in junctions using hybridization expansions [171 09], with 
subsequent work treating the Anderson impurity model [39j and developing expansion in 
the interaction [16]. This first generation of methods was mostly suitable for accessing very 
short times or weakly interacting systems. A much wider range of parameters and timescales 
is accessible to bold-line algorithms [ 7511ST] , which begin from a diagrammatic approxima¬ 
tion containing an infinite subset of diagrams corresponding to a low-order self energy, and 
sum all corrections to it in terms of renormalized skeleton diagrams. These methods can 
be further augmented by reduced dynamics techniques, which give access to essentially any 
timescale in cases where the system exhibits a short memory timescale [43l. HI [82] , More 
recently, these techniques were extended from single-time properties to correlation functions 
in equilibrium [80] and nonequilibrium [83] . 

In this work, we compare our HQME results to bold-line CT-QMC formulated around 
the one-crossing approximation (OCA) [80, 85]. OCA is a strong-coupling approximation. 
It represents an extension of the non-crossing approximation and generally performs well 
near half-filling and outside the Kondo regime. Convergence becomes easier when the OCA 
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is more accurate, but all the CT-QMC data presented here has been converged up to times 
t = 2/r, and can therefore be assumed to be numerically exact, independently of the OCA. 
A detailed technical discussion of the method can be found in Ref. EDI 

While no significant problems occur up to t — 2/r, we note that in general it can be 
difficult to obtain converged CT-QMC data at long times (and in the absence of a short 
memory), since the sign-problem results in an exponential growth of the statistical error 
with time. Bold-line algorithms significantly improve the performance of these algorithms, 
but do not eliminate this problem. “Boldihcation” additionally depends on one’s ability to 
solve the underlying self-consistent diagrammatic approximation (OCA in this case). This 
generally implies an initial step with its own computational and memory demands, both of 
which increase polynomially with the simulated time. Higher-order self-energies reduce the 
sign-problem in the CT-QMC step, but the cost of the diagrammatic approximation often 
becomes prohibitive [83] . 

A particularly simple example, which illustrates how this polynomial scaling can become 
a bottleneck, is when several energy scales which are orders of magnitude apart are present 
in the problem. Unless an efficient multi-scale representation of the data is possible, the 
diagrammatic procedure - which is implemented on a discrete lattice - suffers from the need 
to use very small time steps in the discretization. The effort involved in solving the self- 
consistent equations is then polynomial in the number of time steps, and therefore grows very 
rapidly with simulation time. Importantly, this never occurs with the time-local HQME, 
where the computational effort is inherently linear in time. 


D. Advantages and Drawbacks: When to use which method 

The HQME and CT-QMC schemes are similar in the sense that they are both based 
on a hybridization expansion. The methods differ in the way the expansion is carried out. 
For the HQME approach, we expand the time evolution operator of the reduced density 
matrix. If the expansion converges, one obtains exact results. If not, one obtains only 
approximate results, even at short simulation times. In contrast, the CT-QMC approach 
represents a stochastic sum over all possible trajectories that the system may follow during 
the simulation time. The statistical error or, equivalently, the numerical effort increases in 
the same way as the number of relevant trajectories increases with the simulation time [3lj. 
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The number of relevant trajectories grows exponentially with the simulation time, so QMC 
methods work well for short times, but long-lived correlated dynamics is often out of the 
method’s reach [43]. HQME, in contrast, gives access to long-lived dynamics, because the 
associated numerical effort scales linearly with the simulation time and because all of our 
numerical evidence points out that the quality of the associated expansion is independent of 
the simulation time t. This is, on one hand, evident from Eqs. (Jr]) and, on the other hand, 
explicitly demonstrated in Sec. |III B[ the appendix or, for example, in Ref. 20 We also note 


in passing that because both HQME and CT-QMC are formulated in continous time, they 
are free from discretization (non-zero time step) errors, so that very short times may easily 
be studied. 

We further discuss the numerical effort of the HQME method. Apart from the linear scal¬ 
ing with the simulation time, it depends on the specific problem, in particular the number N 
of distinct superindices j\. The latter is given by the complexity of the correlation functions 
The number of auxiliary operators scales as N K /k\, assuming that the hierarchy ([7]) 


is truncated at the /All tier. The importance criterion (10) further reduces the numerical 
effort to 1V K_1 /R\, cutting out a hypersurface of the total index space [20]. The criterion also 


demonstrates that fewer terms are needed at higher temperatures (cf. the discussion follow¬ 


ing Eq. 10 in Sec. IIB or Ref. 00). In addition, each auxiliary operator involves a number of 
coefficients that is given by the dimension of the Hilbert space of the impurity. In general, 
the size of this space results in an exponential scaling of the numerical effort with the spin 
or orbital degrees of freedom of the impurity. In many cases, however, one is interested in 
single-particle quantities or one can restrict the attention to an active space of considerably 
reduced dimension, possibly enabling a power-law scaling. An explicit demonstration of this 
conjecture will be subject of future research. 

We summarize our discussion regarding the validity and usefulness of the methods in 
Fig. [2j HQME and CT-QMC have common regimes where they give the same result (high 


temperature, short times). This will be demonstrated in Sec. IllB explicitly. There are 
also regimes where only one of the methods can be used in practice. Low temperature 
systems requiring long simulation times cannot be probed by either of the methods, unless 
they also exhibit a short memory timescale, in which case reduced dynamics techniques may 
be applicable. This means that slow dynamics deep in the Kondo regime remain largely 
inaccessible for both methods. The dashed lines in Fig. [2] represent the exponential wall 
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FIG. 2. (Color online) Sketch of the areas in simulation time and temperature where HQME and 
CT-QMC are useful. The dashed lines represent the exponential growth of the numerical effort with 
the simulation time in CT-QMC (black line) and with the number of coefficients in HQME (gray 
line) that increases (at worst exponentially) with the inverse temperature. High temperatures and 
short simulation times are accessible by both methods (white area). Very long simulation times 
are accessible only by HQME (yellow area). The low temperature regime is reserved for CT-QMC 
(blue area). At low temperatures and if long simulation times are required, both CT-QMC and 
HQME cannot be used. For the given problem, the exponential walls are located around the Kondo 


temperature for HQME and time scales ~ 10/T for CT-QMC (cf. Sec. IIIB and Refs. H3l HT , and 


that is hit in CT-QMC with an increasing simulation time and in HQME with the number 
of operators that needs to be taken into account at decreasing temperatures. These walls 
are “soft” in the sense that these boundaries can be pushed by more efficient codes and 
procedures, more powerful computer architectures and merely a larger investment of CPU 
time. They also depend to a large extent on the specific problem. Therefore, we refrain from 
putting specific numbers at this point, but will elaborate on the boundaries specific for the 
Anderson impurity model in the strong coupling regime ( U/(ttT ) > 1) below. 
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E. Observables of interest 


We characterize the nonequilibrium transport properties of an Anderson impurity by its 
magnetization m and the electrical current / that is flowing through the impurity in the 
presence of a bias voltage. The magnetization is given by the diagonal elements of the 
impurity density matrix 

m (t) = Pu(t)~ Putt)’ ( n ) 

where we use the basis {100), | t)> I I), | Tl)}- This basis includes the states of the impurity 
with no electron, 100); a single spin-up electron, | t)i a single spin-down electron, | ),); and 
two electrons, | tl)- 

The electrical current flowing through the impurity is related to the charge flow in and 
out of each lead K: 

= ( 12 ) 

1 k&K 

where —e denotes the charge of an electron. Using the auxiliary operators pp {t), it can be 
written as [19J (52] 

I K {t) =e'PP (Tr imp 

K,(j,p 

III. RESULTS 

A. Formulation of the transport problem 

In the following, we investigate transport and relaxation phenomena of a charge- 
symmetric Anderson impurity where = —U (see Figs, [ljj and 00- We follow 

the time evolution from a product initial state where the impurity is not correlated with 
the electrodes and carries no electron (he. Poo,oo(^ = 0) = 1 while all other elements of the 
reduced density matrix are zero). We focus on the intermediate to strong coupling regime, 
choosing U — 8T (or, equivalently, p ~ 2.5), where T denotes the hybridization strength 
between the impurity and the electrodes at the Fermi level. We have chosen this regime 
because it represents the most challenging regime for the HQME framework and indeed for 
most theoretical treatments (since simple approximations generally work for either T -A 0 or 
T —> oo), while also exhibiting a rich and interesting variety of nonequilibrium phenomena. 




-Tr ; 


imp 


4 pi]*,-,pi*) 


(13) 
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(14) 


We take the tunneling efficiencies (Eq. (J5]) ) to be 


r L /R, ff (e) = ra L /R- 


D- 


[e-Sfx L/R y + D^ 

where we assume Lorentzian-shaped conduction bands in the electrodes with D = 10r. 
Note that the shape of the conduction bands is not crucial for our discussion but beneficial 
for the numerical evaluation of the HQME. The parameters «l/r are either «l = 1 = «r 
corresponding to a symmetric coupling of the impurity to the electrodes or = 1 = 4 «r 
to simulate an asymmetric impurity-electrode coupling. The parameter S is used to control 
whether the conduction bands are shifted with the applied bias voltage (S = 1; cf. Fig. 
[IJ:) or not (S = 0; cf. Fig. [ljo) . While the former situation corresponds to a scenario 
that is typically found, for example, in transport through a single-molecule junction, the 
latter is often used to describe transport through quantum dot structures that are based on 
semiconductor heterostructures. 


Our analysis includes two parts. In the Erst part, Sec. |III B[ we present results for the 
electrical current that is flowing through the impurity in the presence of a bias voltage. 
Thereby, we give a detailed comparison of HQME and CT-QMC results, which, on one 
hand, validates the HQME framework that we introduced in Ref. [T9] (and outlined in 


Sec. IIB) and, on the other hand, allows us to explore the different initial dynamics of the 
electrical current with respect to the two values of the shift parameter S. Second, in Sec. 


(IIC we focus on the dynamics of the dot observable that takes longest to approach its 


steady state value: the magnetization m. We simulate the effect of a magnetic held by 
shifting the spin-up and the spin-down level of the impurity with the held strength h, 


U , 

e t - “7T + h, 


(15) 

(16) 


and study the evolution of m to its field-dependent steady state values. This allows us 
to demonstrate that HQME gives access to long-lived correlated dynamics that, to date, 
had only been accessible at great computational cost using state-to-the-art CT-QMC meth¬ 
ods combined with reduced dynamics [43]. Moreover, we elucidate the origin of the non¬ 
monotonic temperature dependence of the magnetization that was recently reported in Ref. 
[43]. All model parameters are listed in Tab. [TJ where the different parameter sets are labeled 
by the hgure depicting the corresponding results. 
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TABLE I. Parameters for the quantum dot devices that are investigated in this article. All energy 
values are given in units of the hybridization strength T. 


B. Time-dependent electrical current: Comparison of HQME and CT-QMC 


The primary goal of this section is to provide the first direct comparison of results that 
have been obtained by the HQME and CT-QMC approach. We thus validate the HQME 
scheme with respect to an established method and, at the same time, demonstrate explicitly 
that our truncation scheme is consistent and gives numerically exact results once convergence 


is achieved. As the importance criterion (10) suggests, the HQME expansion works best for 
high temperatures. Accordingly, we start our comparison at a relatively high temperature 
in the electrodes, /3 = = 0.2/T, and continue with an intermediate, f3 = 1/T, and 

a low temperature, f3 = 5/T. The latter is close to the Hondo temperature of our setup, 

/^Kondo ~ 15/T 


We begin with the case S = 0, corresponding to fixed bands in the electrodes. Fig. 
[3] represents the symmetrized current / = (Jl — /r)/ 2 flowing through the impurity as a 
function of time at /3 = 0.2/T. The different lines correspond to bias voltages e<f> = T, 3T, 
... 19T. The thick blue and thin orange lines depict results that have been obtained using 
the HQME and CT-QMC scheme, respectively. The overlap between matching pairs of lines 
demonstrates the agreement between HQME and CT-QMC in this range of temperatures. 

The dynamics seen in Fig. [3] can be described and understood as follows. Initially, for 
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t < 0.2/r, the current increases almost linearly in time with a slope that is increasing linearly 
with the applied bias voltage. After this initial increase, the current level saturates rapidly 
to a stationary state. The short-time behavior can be understood quantitatively from the 
relation 

^n-(o) = 2eJ2 (c£,„(ox<i„4> - cy„(o)<4<f„», ( 17 ) 

a 

which follows directly from the operator equations of motion and the choice of a product 

initial state. Ck,<t is defined in Eq. ([8]). For an initially unoccupied impurity, the slope of 

the symmetrized current is therefore given by 
d t'°° d'; 

Tt I{t) = 2e£ J ^ • (18) 

For large band width D, the energy dependence of the hybridization strengths Tr )(T (u;) can 
be neglected. The slope of the current is then solely determined by the difference between 
the Fermi functions. The latter is proportional to the applied bias voltage <f>. Note that the 
initial dynamics cannot be linked to a single time scale here but is influenced by the position 
of the energy levels, the band-width D and the temperature T. 

At lower temperatures, both schemes require a larger computational effort in order to 
reach the same level of precision as compared to higher temperatures. This is demonstated 
in Figs. [4] and [5j which show the time-dependent current of our setup at lower temperatures 
P = 1/T and /3 = 5/T, respectively. The data has been obtained with a similar numerical 
effort as that shown in Fig. [3j One observes that both schemes agree very well, but small 
deviations, which are consistent with the applied accuracy, begin to occur. 

As the temperature T decreases, coherent processes become more important and give 
rise to oscillations of the current level (see Figs. [4] and [5]). The period of these oscillations 
is given by the dynamical phases of the system, in particular the difference between the 
energy levels of the impurity and the chemical potentials in the electrodes. Accordingly, the 
dynamical oscillations of the current level show a bias dependence, which is clearly visible 
in our data. This behavior is known as current ringing and has been outlined before in 
a slightly different context, namely as a response to bias voltage pulses and/or quenches 

[STIES]. 

Another bias dependence appears in the stationary values of the current level. For high 
temperatures, thermal broadening leads to an almost linear increase of the stationary cur¬ 
rent, at least in the range of bias voltages considered here. This is evident from the almost 
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FIG. 3. (Color online) Symmetrized current / = (if, — /r)/2 flowing through the impurity at 
k^T = 5r as a function of time t for a sequence of equally spaced bias voltages <f> = T, 3r,.., 19r 
where the conduction bands are not shifted with the bias voltage (S = 0). The model parameters 
used to obtain this data are summarized in Tab. [T] After a linear increase, the current saturates to a 
stationary value on a voltage-independent time scale 0.2/T. The HQME (blue lines) and CT-QMC 
methods (orange lines) give identical results to within the numerical resolution of the data. 

equidistant values in Fig. [3j The stationary values seen in Figs. [| and [5] are clearly non- 
equidistant. This indicates a strong non-Olnnic saturation of the current level with increasing 
bias voltage, which originates from the restricted number of conductance channels through 
the impurity. 

We conclude at this point that the agreement between HQME and CT-QMC results is 
very good in the parameter ranges we have studied. We corroborated this statement for a 
number of other setups, where we changed the position of the energy levels, ^ —U/ 2, 
introduced a level splitting / magnetic held, G — ^ 0, a different shift of the conduction 

bands, S — 1 (see below), and different band widths (data not shown). Changing the 
electron-electron interaction strength, we also observed that the results converge faster for 
lower values of U (which is consistent with the fact that, for U — 0, the hierarchy ([7]) 
terminates at the second tier mm- At lower temperatures, f3 > /^Rondo ~ 15/T, we were 
not able to converge the HQME expansion to a satisfactory level. 

Next, we consider electrodes in which the conduction bands are shifted with the applied 
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FIG. 4. (Color online) Symmetrized current / = (/l — /r)/ 2 flowing through the impurity at 
k^T = T as a function of time t for a sequence of equally spaced bias voltages = T, 3r,19r 
where the conduction bands are not shifted with the bias voltage (S = 0). The model parameters 
used to obtain this data are summarized in Tab. [I] After a linear increase, the current level slightly 
oscillates before it reaches its stationary value. The corresponding time scales are given by the 
dynamical phases of the problem and the inverse temperature 1 /(/cbT) = 1/T, respectively. 


bias voltage (S = 1). The time-dependent current of such a system is shown in Fig. [6j 
It should be compared and contrasted with the data shown in Fig. |4j Two qualitative 
differences are apparent: the first is that the slope of the current vanishes at t = 0. This can 


be easily understood from Eq. (18), as the integral on the right-hand side vanishes for shifted 
conduction bands. The latter is no longer true for times t > 0, where, initially, an increase 


of the current ~ t 3 is inherited from a change of the populations ~ t 2 


After 


this initial phase, the behavior of the unshifted bands is recovered. The other difference 
is the reduced current level for e<F = 19T, which falls below the line for e<F = 11T at 
times t > 0.8/r. This negative differential resistance originates from both the shift of the 
conduction bands and their finite band width D. This behavior is also well known, for 
example, from transport through resonant tunneling diodes 
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FIG. 5. (Color online) Symmetrized current / = (Jl — /r)/ 2 flowing through the impurity at 
k^T = T/5 as a function of time t for a sequence of equally spaced bias voltages <J> = T, 3r,19r 
where the conduction bands are not shifted with the bias voltage (S' = 0). The model parameters 
used to obtain this data are summarized in Tab. [TJ The oscillations of the current level that appear 
right after the initial linear increase become more pronounced at lower temperatures. The corre¬ 
sponding time scales are given by the dynamical phases of the problem and the inverse temperature 
l/(kftT) = r/5, respectively. Deviations between the HQME (blue lines) and CT-QMC results 
(orange lines) are consistent. 


C. Evolution of the magnetization for a symmetric and an asymmetric coupling 
to the electrodes 


The HQME method is particularly promising because it allows us to calculate the exact 
time evolution of a correlated many-body system with a numerical effort that scales linearly 
with the simulation time. In order to demonstrate this, we discuss the magnetization, as 


introduced in Eqs. (15). Other observables, like populations or the current, approach their 


steady state values much faster and are, therefore, less suitable for the present purpose. 
In addition, we consider an asymmetric coupling to the electrodes. As we will see, the 
corresponding dynamics shows a richer set of behaviors than in the symmetrically coupled 
case, and occurs on significantly longer time scales. 

We start our analysis with a symmetrically coupled impurity in a strong magnetic field 
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FIG. 6. (Color online) Symmetrized current / = (if, — /r)/ 2 flowing through the impurity at 
k-^T = T as a function of time t for a sequence of equally spaced bias voltages = T, 3r,lir 
and = 19r where the conduction bands are shifted with the bias voltage (S = 1). The model 
parameters used to obtain this data are summarized in Tab.[IJ In contrast to Figs. !-! the current 
level is initially not linearly increasing with the time t when the conduction bands are shifted with 
the applied bias voltage (S = 1). 


h = 2T (the behavior at other held strengths is similar - data not shown). The corresponding 
magnetization of the impurity is depicted in Fig. [7] as a function of both bias voltage <f> and 
inverse temperature in the electrodes f3. We show a 3 x 3 array of plots, where the top 
row represents the magnetization at short time scales, t = 0.5/T, and the middle row at 
intermediate time scales, t = 5/T. The bottom row depicts the steady state values. In 
the latter, we also give the time fo .999 at which the impurity reached 99.9% of its final 
magnetization (note that this time scale is longest for low temperatures and bias voltages). 
The different columns are obtained using different levels of approximation. The left column 
is computed using the full HQME approach. The second and third column are obtained 
by truncating the hierarchy (J7]) at the second and the first tier. This corresponds to a 
hybridization expansion to 0((T / (ksT)) 2 ) and 0(T/(k^T)), respectively (see Sec. IIB). 
The different levels of approximation facilitate a discussion of the relevant processes and 
mechanisms, for example, where and when processes of higher order are important (see 
below). 
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FIG. 7. (Color online) Magnetization of an Anderson impurity that is symmetrically coupled to 
the electrodes as a function of temperature and bias voltage (5 = 0), and for three different times 
(top row: t = 0.5/r; middle row: t = 5/T; bottom row: steady state). The left column depicts the 
full HQME result. The middle and the right column has been obtained by truncating the hierarchy 
([7j) at the second and the first tier, respectively. The magnetization along the red dashed line in 
the lower left plot is also depicted in Fig. |8j 


The data of Fig. [7] can be understood as follows. As we start from a state with zero 
magnetization, the magnetization at short time scales is almost an order of magnitude 
smaller than the final steady state values. The maximum magnetization is obtained at small 
bias voltages and temperatures. At higher temperatures and voltages, where the impurity 
exchanges particles with the electrodes at a wider range of energies, the magnetization 
becomes quenched. We would like to emphasize that the system reaches its steady state not 


21 














































































































FIG. 8. (Color online) Magnetization of an Anderson impurity that is symmetrically coupled to 
the electrodes as a function of temperature in the steady state and applied bias voltages = ±13r, 
±15r and ±17r (S = 0). Note that the kinks originate from a simple linear interpolation between 
the data points. 

before times t > 15/T. 

The exact result is very similar to the one that is obtained with a second order truncation 
of the hierarchy ([7]). A tendency towards a higher magnetization is visible if the number 
of exchange processes that is taken into account in our calculations is reduced. Truncation 
at the first tier enhances the effect, but also results in a qualitatively different structure 
of the magnetization at intermediate time scales (right plot of the middle row). Here, the 
magnetization shows a peak at positive and negative voltages, while the exact and second 
order result exhibit only a single peak that is centered around zero bias. 

The splitting of the peak magnetization can be understood with the bias dependence 
of resonant exchange processes with the electrodes. At zero bias, the difference between 
the chemical potentials in the electrodes and the single-particle levels of the impurity (see 
Fig. [I|d) is largest and resonant processes are strongly suppressed. This suppression is less 
pronounced at non-zero voltages such that the steady state magnetization can develop on 
shorter time scales. Accordingly, this behavior shows only a weak temperature dependence, 
resulting in an almost horizontal splitting of the peak that is seen in Fig. [TJ This splitting 
is not seen at short times scales (right plot of the top row), because the initial state is not 
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decaying exponentially at short times. It rather shows a power-law decay, ~ t 2 , which is 
well known from an analysis of similar systems in terms of Born-Markov theory [SO, [80J 100 J. 

Another intriguing effect occurs at higher bias voltages (|e<h| > 10r). Here, the magneti¬ 
zation shows a non-monotonic behavior with respect to temperature: it becomes stabilized 
by increasing temperature before decreasing again at temperatures f3 < 0.2/r (follow, e.g., 
the red dashed line in Fig. [7] from /3 = 1/T to /3 = 0.1/r). This non-monotonic behavior is 
explicitly depicted in Fig. [8j where the magnetization m is shown, for example, as a function 
of the inverse temperature f3 and fixed bias voltages <f> = ±13r, ±15r and ±17r. This non¬ 
linear dependence of the magnetization on temperature was discovered only recently (see 
Ref. |43j). It is most pronounced in the steady state and for a truncation of the hierarchy 
0 at a lower tier. The latter points towards a physical interpretation of the effect, suggest¬ 
ing that it originates from the broadening of the peak magnetization around zero bias with 
increasing temperature and the quenching of the magnetization at very high temperatures. 

The last scenario we discuss is an asymmetric coupling of the impurity to the electrodes. 
The corresponding magnetization is shown in Fig. [9j It can be directly compared to the 
magnetization of the symmetrically coupled impurity that is depicted in Fig. [0 It can be seen 
that, on short time scales (top row), the initial peak of the magnetization is shifted towards 
negative voltages. This is because the coupling to the right electrode is weaker and we start 
with an initially unoccupied system. The initial population of the impurity is therefore 
dominated by exchange processes with respect to the left electrode. The corresponding 
dynamics occurs on shorter time scales for negative voltages, because the chemical potential 
of the left electrode /il is then closer to the single-particle levels e-pq. 

Another qualitative difference with respect to the symmetric case occurs on intermediate 
time scales. Here, the magnetization is peaked at non-zero values of the bias voltage <f> 
even when higher-order processes are taken into account. This behavior was also seen 
in the symmetric case but only if the hierarchy 0 is truncated at the first tier, that is 
by disregarding higher-order processes. These processes become quenched by the weaker 
coupling to the right electrode, while resonant processes with respect to the left electrode 
are not. The situation is therefore similar to the symmetrically coupled case without higher- 
order processes. The magnetization of the impurity evolves to very similar steady state 
values, but on even longer times scales (t > 30/T). Minor differences occur due to the 
weaker hybridization with the right electrode [i.e. a less pronounced broadening of the 
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FIG. 9. (Color online) Magnetization of an Anderson impurity asymmetrically coupled to the 
electrodes as a function temperature and bias voltage (S = 0) and for three different times (top 
row: t = 0.5/r; middle row: t = 5/T; bottom row: steady state). The left column depicts the full 
HQME result. The middle and the right column has been obtained by truncating the hierarchy 
|7| at the second and the first tier, respectively. 


energy levels). 

We close this section by a discussion on the generality of our hirelings. We observed, 
for example, a very similar behavior of the magnetization dynamics for different choices of 
the voltage division factor (/tl ^ — /ir for $ ^ 0, data not shown). We also calculated 
the magnetization dynamics starting from different initial states. If, for example, the ini¬ 
tial magnetization points in the direction of the magnetic held, the impurity magnetization 
shows a very similar behavior as compared to the one we discussed for a symmetric cou- 
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pling to the electrodes. Similar structures as for an asymmetric coupling to the electrodes 
appear, if the initial magnetization points opposite to the applied magnetic field. These 
signatures, however, decay on much shorter time scales due to the stronger coupling to the 
right electrode. 


IV. CONCLUSION 


In this work, we give the first direct comparison of the hierarchical quantum master 
equation method m 1201 EQH52] and the diagrammatic, continuous-time quantum Monte 
Carlo approach [371 ESI 031 mi SHI m]. To this end, we have studied the nonequilibrium 
transport properties of an Anderson impurity that is coupled to two electrodes with different 
chemical potentials. This transport problem represents a well established and fairly well 
understood test case. We discussed the main characteristics of the two numerically exact 


methods (cf. Sec. IID). They are distinguished by the range of parameters where exact 


results can be obtained in practical calculations. CT-QMC gives access to the short- and 
intermediate-time dynamics (< 10 units of the inverse hybridization strength) but in general 
fails to describe long-time dynamics, for example in the presence of Kondo correlations, 
because the numerical effort scales exponentially with the simulation time [32]• In contrast, 
the hierarchical master equation method scales linearly with the simulation time and is, 
therefore, not limited in terms of the accessible time scale. It represents a hybridization 
expansion, which can be carried out to sufficiently high order if the temperature in the 
electrodes is not too low. For the present problem, we were able to obtain converged results 
only for temperatures above the Kondo temperature. Our findings provide a wealth of 
numerically exact benchmark data certain to be useful to future method developers. 

We have also elucidated interesting physical phenomena for the range of parameters, 
where numerically exact results have been accessible with a reasonable numerical effort. 
We investigated the short-time dynamics of the (symmetrized) current flowing through the 
impurity in the presence of a bias voltage, starting from a product initial state where the 
impurity is not populated by electrons. While a linear increase of the current level is found 
for situations where the conduction bands are not shifted with the applied bias voltage 
(corresponding to realizations of quantum dots with semiconductor heterostructures), a 
qualitatively different behavior emerges when the electrodes are not charged when applying 
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a bias voltage. At low temperatures, oscillations of the current level, current ringing [87], due 
to dynamical phases appear. In addition to the current, we also studied the magnetization 
dynamics in the presence of a magnetic field. We confirmed the recently reported non¬ 
monotonic temperature dependence of the steady state magnetization [33] and traced it 
back to competing broadening effects of the impurity levels. In addition, we found complex 
structures on intermediate yet long time scales (i.e., in the cases studied, about 10 units 
of the inverse hybridization strength) in the presence of an asymmetric coupling to the 
electrodes (cf. Fig. [9]). 

Our comparative study is a first step towards practical guidelines in choosing the right 
solver for a given impurity problem. Since the presented methods cannot cover the full 
spectrum of problems, it would be interesting to compare them to other exact schemes, 
including, for example, numerical renormalization group, density matrix renormalization 
group, multi-layer multi-configurational time-dependent Hartree, or other quantum Monte 
Carlo schemes. A primary goal is to identify regions of parameter space where methods 
overlap and, of course, regions which cannot be covered satisfactorily by any of the available 
methods. Further activities in this direction are planned, and we would like to encourage 
other researchers to participate in these efforts. 
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Appendix A: Convergence properties of HQME 

The convergence properties of the HQME approach strongly depend on the importance 
criterion that is used to truncate the hierarchy of equations of motion ([7]) . In this appendix, 
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threshold value 

10” 1 

10~ 2 

10~ 3 

10 -4 

10 -5 

10 -6 

# of AOs 

269 

447 

1158 

3009 

8317 

20912 

max. tier level 

2 

3 

3 

4 

5 

5 


TABLE II. Number of auxiliary operators for different threshold values of the importance criterion 


( 10 ). 


we demonstrate the convergence behavior of our HQME scheme explicitly. To this end, we 
reconsider the time-dependent (converged) current shown in Fig. |4j We replot the result for 
$ = 5T in Fig. 10 1 . It corresponds to the graph with the threshold value A t h = 10~ 6 . In 
addition, we plot results that were obtained for higher threshold values A t h- In Tab. [TT} we 
list the number of auxiliary operators (AO) that were taken into account and the maximum 
tier level. The results are considered to be converged once the threshold value is below 10~ 5 . 
The corresponding number of AOs is ~ 10 4 . The respective maximum tier level is 5. 


In Sec. Ill B we have shown that our (converged) results coincide with the ones obtained 
from CT-QMC. This demonstrates the validity and usefulness of our truncation scheme. 
At this point, we would like to give an additional proof of this statement by showing the 
convergence of our scheme towards the solution of an analytically solvable case. Fig. [Id} ) has 
been obtained with the same parameters as Fig. [To) ), except that we ’turned off’ electron- 
electron interactions, i.e. we used U = 0. In this limit, the exact result is known and can be 
obtained, for example, by truncating the hierarchy ([?]) at the second tier (without applying 


the importance criterion (10)) [52j 172] . It can be seen that our results converge to the exact 
result and that convergence is faster than in the interacting case. 
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FIG. 10. (Color online) Convergence analysis of the symmetrized current I = (Jl — /r)/ 2 that is 
flowing through the impurity at k^T = T and <h = 5r (S = 0) as a function of time t. Panel (a) 
shows the convergence behavior for U = 8r. The behavior for U = 0 is shown in Panel (b). The 
model parameters used to obtain this data are summarized in Tab. [T| 
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